River Hydrodynamics Simulation

Hydraulic infrastructures management and operation are challenging issues. Considering climate changes, risks and uncertainties, it is increasingly important to assure environmental, economic and social sustainability. The use of adequate and efficient tools can contribute with solutions supported on hydrodynamic models, used to simulate the flow of water along river channels and floodplains.

On the other hand, Web-based tools can provide a suitable environment to illustrate the hydrodynamics behavior of a given river under different conditions and to evaluate various scenarios and management strategies. Moreover, this approach can also be considered as an academic resource for learning purposes.

This work presents a Web interface, supported by a Jupyter Notebook, which provides functionalities to configure a hydrodynamic river model and simulate its behavior. The model setup is run through the MIKE 1D numerical simulation engine behind MIKE HYDRO River commercial software.

This approach considers the Mondego River as a case study and allows the use of different datasets for different meteorological conditions, although here, for demonstration purposes, a particular dataset is used. This tool also provides access to real data, such as water levels and streamflows, collected by geosensors located in the considered area.

Throughout this notebook, the successive tasks are explained. Let's get started by loading some Python libraries.

In [1]:
%load_ext autoreload
%autoreload 2
In [2]:
from pathlib import Path
import pandas as pd

import geoviews as gv
import geopandas as gpd
import cartopy
from cartopy import crs as ccrs

import holoviews as hv
import hvplot
import hvplot.pandas
import panel as pn
import param

import src.riverhydrosim.mike1Dsimulation as sim
import src.riverhydrosim.utils as utils
In [3]:
# pandas options
pd.set_option('display.max_rows', 10)
In [4]:
pn.extension()
hv.extension('bokeh')

Define the folders where are the datasets and the river model to run:

In [5]:
MODELS_DIR = "./models/Mondego_Simulacao_CheiaJan2016"
DATA_DIR = "./data"

1. Case study - the Portuguese Mondego River

The Mondego River is the fifth biggest Portuguese river, being the largest river with both headwaters and mouth on Portuguese territory, located in central region of Portugal.

The case study consists of the stretch of Mondego River, between the Aguieira-Raiva dams and the dam bridge of Coimbra City, with an approximately 36.5 km long. The case study considers the two main tributaries, Alva River and Ceira River. The streamflow of the Alva River is controlled by the Fronhas Dam, but the Ceira River streamflow is uncontrolled.

The Mondego River widens substantially around Coimbra City making it a flood-prone area. Since the riparian zones of Coimbra are often affected by river floods at more intense rainfall events, the study of the river hydrodynamics is fundamental for flood forecasting and risk analysis and will help decision-makers to be better prepared for floods.

Some images of past flood events (January 2016) in Coimbra city:

Drawing
Santa Clara Convent
Drawing
Green Park

The location of the river, the key cross sections and the sensors is shown on the map.

In [6]:
river_shp = Path("./shapefiles/Mondego_opt_Branches.shp")
cross_sections_shp = Path("./shapefiles/Mondego_opt_CrossSections.shp")
inflow_points_shp = Path("./shapefiles/inflow_points.shp")
sensors_shp = Path("./shapefiles/sensors.shp")

opts_1 = dict(height=500, width=700, fontsize={"legend": 8, "title": 12}, active_tools=["box_zoom", "pan"])
opts_2 = dict(show_legend=True, apply_ranges=True, tools=["hover"])
# basemap = gv.tile_sources.OSM
basemap = gv.tile_sources.EsriImagery()

# river shape
r = gpd.read_file(str(river_shp), encoding="ISO-8859-1")
river = gv.Path(r, crs=ccrs.GOOGLE_MERCATOR, vdims=[("BR_BrName", "branch_name")], 
                label="Mondego river").opts(color="blue", line_width=1.5, **opts_2)

# key cross section shapes 
cs = gpd.read_file(str(cross_sections_shp), encoding="ISO-8859-1")
key_cs = cs.iloc[[15,16,18,20]]
cross_sections = gv.Path(key_cs, crs=ccrs.GOOGLE_MERCATOR, 
                         vdims=["CS_ID", ("CS_Ch", "CS_chainage")],
                         label="Cross sections").opts(color="red", line_width=2, **opts_2)

# boundary (inflow) point shapes
inflows = gpd.read_file(str(inflow_points_shp), encoding="ISO-8859-1")
inflow_points = gv.Points(inflows, crs=ccrs.GOOGLE_MERCATOR, 
                          vdims=["Name", ("UDC_Ch", "BR_chainage"), ("UDC_BrName", "branch_name")],
                         label="Boundaries (inflows)").opts(color="cyan", size=7, **opts_2, legend_position="top_left")

# sensor shapes
s = gpd.read_file(str(sensors_shp), encoding="ISO-8859-1")
sensors = gv.Points(s, crs=ccrs.GOOGLE_MERCATOR, 
                    vdims=[("Code_3", "SE_ID"), ("Name", "name"), ("Code", "SNIRH_code")], 
                    label="Sensors").opts(color="yellow", size=7, **opts_2, legend_position="top_left")

p = basemap.opts(**opts_1) * river * cross_sections * inflow_points * sensors
p.opts(padding=2.0)
Out[6]:

2. Loading and viewing geosensor observations

There are two options for loading data from geosensors:

  1. Remotely, from SNIRH - Sistema Nacional de Informação de Recursos Hídricos database;
  2. Or locally, from a previously downloaded CSV file.

Let's load some data observations (water levels), for instance between 2016-01-01 and 2016-01-12, when Coimbra city suffered a severe flood inundation (see images above).

In [7]:
tmin="2016-01-01"
tmax="2016-01-12"

# REMOTE access
args = ("remote", "", tmin, tmax)

# LOCAL access
# args = ("local", Path(DATA_DIR) / "processed/water_level_observ_20160101_20160112_.csv", tmin, tmax)

def load_water_level_observations(source="local", file="", tmin="", tmax=""):
    tmin = pd.to_datetime(tmin)
    tmin = tmin.strftime(format="%d-%m-%Y")
    tmax = pd.to_datetime(tmax)
    tmax = tmax.strftime(format="%d-%m-%Y")
    
    try:
        if source == "remote":
            info = {"code": ["12G/01A", "12G/04H"],
                "name_abbrev": ["APC", "PSC"],
                "sites": [1627743374, 1627759222],
                "pars": [1629599726, 1843],
                "ref_level": [0.0, 15.45]
               }

            url_template = ("https://snirh.apambiente.pt/snirh/_dadosbase/site/paraCSV/"
                   "dados_csv.php?sites={}&pars={}&tmin={}&tmax={}&formato=csv")
            dfs = []
            for i in range(len(info["name_abbrev"])):
                url = url_template.format(info["sites"][i], info["pars"][i], tmin, tmax)
                df_obs = pd.read_csv(url, skiprows=4, encoding="latin-1")
                if df_obs.empty:
                    print("Unable to download water level observations from sensor {}.".format())
                else:
                    # print(url)
                    df_obs = pd.read_csv(url, skiprows=4, names=["datetime", "water_level"], 
                                         usecols=[0,1], encoding="latin-1")
                    df_obs = df_obs[:-1] # last row is garbage
                    # format datetime column
                    df_obs["datetime"] = pd.to_datetime(df_obs["datetime"], format='%d/%m/%Y %H:%M')
                    df_obs["water_level"] = df_obs["water_level"] + info["ref_level"][i]
                    df_obs["name_abbrev"] = info["name_abbrev"][i]
                    df_obs["code"] = info["code"][i]
                    dfs.append(df_obs) 

            if dfs: # list not empty
                return pd.concat(dfs, ignore_index=True) # append all dataframes in dfs list
            else:
                return None #pd.DataFrame() # empty dataframe

        elif source == "local":
            return pd.read_csv(file, usecols=["datetime", "water_level", "name_abbrev", "code"], encoding="latin-1", 
                               parse_dates=["datetime"], infer_datetime_format=True)
    except Exception as e:
        print(e)

df_obs = load_water_level_observations(*args)
df_obs.rename(columns={"water_level": "water_level (m)", "name_abbrev": "sensor_ID", "code": "SNIRH_code"})
Out[7]:
datetime water_level (m) sensor_ID SNIRH_code
0 2016-01-01 00:00:00 17.93 APC 12G/01A
1 2016-01-01 01:00:00 17.91 APC 12G/01A
2 2016-01-01 02:00:00 17.90 APC 12G/01A
3 2016-01-01 03:00:00 17.89 APC 12G/01A
4 2016-01-01 04:00:00 17.89 APC 12G/01A
... ... ... ... ...
571 2016-01-12 19:00:00 18.65 PSC 12G/04H
572 2016-01-12 20:00:00 18.64 PSC 12G/04H
573 2016-01-12 21:00:00 18.63 PSC 12G/04H
574 2016-01-12 22:00:00 18.62 PSC 12G/04H
575 2016-01-12 23:00:00 18.62 PSC 12G/04H

576 rows × 4 columns

Plot geosensor observations between 2016-01-01 and 2016-01-12

In [8]:
opts = dict(width=700, height=400, fontsize={"legend": 8, "title": 10}, legend_position="top_left")
plot_sensors_data = df_obs.hvplot(x="datetime", y="water_level", ylabel="water level (m)", 
                  title="Geosensor observations: water level", by="name_abbrev",
                  line_width=1.5, color=["green", "orange"]).opts(**opts)
plot_sensors_data
Out[8]:

3. Loading and viewing the system inflows

System inflows are located at:

  • Aguieira-Raiva dams
  • Alva River mouth
  • Ceira River mouth (+ intermediate watershed)

The inflow time series studied in detail, begin at 2016-01-01 00:00:00 and end at 2016-01-12 23:00:00. They can be loaded from a spreadsheet or a CSV file.

In [9]:
# Data file
data_file = Path(DATA_DIR) / "raw" / "data_mondego_2016.xlsx"
# data_file = Path(DATA_DIR) / "raw" / "data_mondego_2016.csv"
cols = dict(zip(["data_hora", "Q_efluente_RV", "Q_efluente_FR", "Q_ceira_bacia_interm"], 
                ["datetime", "Qin_Aguieira_Raiva", "Qin_Alva", "Qin_Ceira"]))

# inflows dataframe
df_inflows = pd.read_excel(data_file, sheet_name="cheia_jan2016", usecols=cols.keys(), 
                           index_col="data_hora").loc["2016-01-01 00:00:00":"2016-01-12 23:00:00"]
df_inflows.index.rename(cols["data_hora"], inplace=True)
df_inflows_ = df_inflows.rename(columns={k: v + " (m3/s)" for (k, v) in cols.items()}).round(2)
df_inflows.rename(columns=cols, inplace=True)
df_inflows_
Out[9]:
Qin_Aguieira_Raiva (m3/s) Qin_Alva (m3/s) Qin_Ceira (m3/s)
datetime
2016-01-01 00:00:00 0.00 2.00 7.92
2016-01-01 01:00:00 0.00 2.00 10.26
2016-01-01 02:00:00 0.00 2.00 11.19
2016-01-01 03:00:00 0.00 2.00 8.51
2016-01-01 04:00:00 0.00 2.00 16.67
... ... ... ...
2016-01-12 19:00:00 735.93 53.62 186.89
2016-01-12 20:00:00 733.80 53.60 181.59
2016-01-12 21:00:00 733.66 53.60 179.67
2016-01-12 22:00:00 732.72 53.60 184.76
2016-01-12 23:00:00 731.93 53.58 174.96

288 rows × 3 columns

Calculate some basic statistics of the inflows

In [10]:
df_inflows.describe()
Out[10]:
Qin_Aguieira_Raiva Qin_Alva Qin_Ceira
count 288.000000 288.000000 288.000000
mean 275.395104 17.943472 117.150333
std 278.065287 22.639282 92.248113
min 0.000000 2.000000 1.968570
25% 0.000000 2.000000 22.786609
50% 159.690000 2.000000 102.407000
75% 431.887500 51.552500 180.149000
max 1108.390000 53.640000 425.098000
In [11]:
opts = dict(width=700, height=400, fontsize={"legend": 8, "title": 10}, legend_position="top_left")
plot_inflows = df_inflows.hvplot(title="Inflows", xlabel="datetime", ylabel="discharge (m3/s)").opts(**opts)
plot_inflows
Out[11]:

4. Convert the inflows in CSV/spreadsheet format to dfs0 format

In order to run the MIKE 1D numerical simulation engine, the time series inflows must be in a proprietary format - dfs0. So, CSV/spreadsheet files must be converted first.

In [12]:
output_folder = Path(DATA_DIR) / "processed"
output_files = [col + ".dfs0" for col in df_inflows.columns]
for i, f in enumerate(output_files):
    dfs0_file = output_folder / f 
    ds_discharge = df_inflows[output_files[i].split(".")[0]]
    utils.convert_discharge_ts_to_dfs0_2(ds_discharge, dfs0_file)

Let's verify if the files were actually created:

In [13]:
! ls ./data/processed/Qin*
./data/processed/Qin_Aguieira_Raiva.dfs0
./data/processed/Qin_Alva.dfs0
./data/processed/Qin_Ceira.dfs0

5. Configure hydrodynamic river model parameters, simulate and visualize the results

A 1D simplified hydrodynamic river model was developed using the software MIKE HYDRO River (not shown in this notebook). This software has a relatively complex graphical interface. This type of model requires the configuration of many parameters and the schematic definition of the fluvial network based on a map.

After the model setup file has been created, it can be run through the MIKE 1D numerical simulation engine behind MIKE HYDRO River.

MIKE 1D is a 1-dimensional numerical simulation engine that simulates water flow in a network (e.g., river, water distribution, storm water drainage and sewer collection networks). The MIKE 1D API provides a high level of interaction and control of the MIKE 1D engine and makes possible to programmatically create or modify a setup and interact with the engine during simulation from user defined code. The same functionality available in the setup editors can be used at runtime.

Here we only want to adjust some of the parameters and analyze the effect on the results - water levels and discharges - in some of the cross sections of the river, those where flooding usually occurs. The key cross sections, located in Coimbra City flood-prone zone, are:

  • Clube naútico (CN)
  • Choupalinho (CH)
  • Ponte Santa Clara (PSC)
  • Açude-ponte Coimbra (APC)

The graphs are interactive and synchronized with each other. The user can exploit different sets of parameters values and quickly inspect the results. The notebook code can be extended to perform other relevant data analysis, such as statistical ones.

In [14]:
class Mike1DSimulation_parameters(param.Parameterized):
    # parameterized inputs
    dt_min = pd.datetime(2016, 1, 1, 0, 0, 0)
    dt_max = pd.datetime(2016, 1, 12, 23, 0, 0)
    simulation_start = param.Date(default=pd.datetime(2016, 1, 5, 0, 0, 0), 
                                  bounds=(dt_min, dt_max),
                                  doc="Simulation start date/time.") 
    simulation_end = param.Date(default=pd.datetime(2016, 1, 12, 23, 0, 0), 
                                bounds=(dt_min, dt_max),
                                doc="Simulation end date/time.")
    time_step = param.Integer(default=5, bounds=(1,30), step=1)
    time_step_unit = param.ObjectSelector(default="minutes", objects=["minutes", "seconds"])
    storing_frequency = param.Integer(default=60, bounds=(1,60), step=5)
    storing_frequency_unit = param.ObjectSelector(default="minutes", objects=["minutes", "hours"])
    grid_spacing = param.Integer(default=1000, bounds=(50,2000), step=50, doc="Grid spacing in meters.")
    delta = param.Number(default=0.5, bounds=(0.5, 1.0), step=0.05, doc="Delta value.")
#     button = param.Action(label="Run")

    # other inputs
    hydro_model_setup_file = str((Path(MODELS_DIR) / "Mondego_opt.mhydro").resolve())
    verbose = False
    chainages = ('clube_nautico', 'choupalinho', 'pte_sta_clara', 'acude_pte_cbra')
    dfs0_aguieira =  str((Path(DATA_DIR) / "processed" / "Qin_Aguieira_Raiva.dfs0").resolve())
    dfs0_foz_alva = str((Path(DATA_DIR) / "processed" / "Qin_Alva.dfs0").resolve())
    dfs0_foz_ceira_baciainterm = str((Path(DATA_DIR) / "processed" / "Qin_Ceira.dfs0").resolve())
    
    

    @param.depends("time_step_unit", watch=True)
    def _update_time_step_bounds(self):
        if self.time_step_unit == "minutes":
            self.param.time_step.bounds = (1, 30)
            self.param.set_param(time_step=5)
        elif self.time_step_unit == "seconds":
            self.param.time_step.bounds = (1, 60)
            self.param.set_param(time_step=30)

    @param.depends("storing_frequency_unit", watch=True)
    def _update_storing_freq_bounds(self):
        if self.storing_frequency_unit == "minutes":
            self.param.storing_frequency.bounds = (1, 60)
            self.param.storing_frequency.step = 5
            self.param.set_param(storing_frequency=60) 
        elif self.storing_frequency_unit == "hours":
            self.param.storing_frequency.bounds = (1, 6)
            self.param.storing_frequency.step = 1
            self.param.set_param(storing_frequency=1)
            
    # =================================================================================================================            

    # run MIKE 1D simulation and plot results (water level and discharges in key cross sections)
    def run_and_plot_results(self):
        try:
            # run MIKE 1D simulation and plot results (water level and discharges in key cross sections)
            df = sim.run_simulation(hydro_model_setup_file=self.hydro_model_setup_file,
                                    chainages=self.chainages,
                                    simulation_start=str(self.simulation_start), 
                                    simulation_end=str(self.simulation_end),
                                    delta=self.delta, grid_spacing=self.grid_spacing,
                                    time_step=(self.time_step, self.time_step_unit),
                                    time_span=(self.storing_frequency, self.storing_frequency_unit),
                                    dfs0_aguieira=self.dfs0_aguieira,
                                    dfs0_foz_alva=self.dfs0_foz_alva,
                                    dfs0_foz_ceira_baciainterm=self.dfs0_foz_ceira_baciainterm,
                                    verbose=self.verbose)
            
            self.results = df

            # build water level and discharge charts
            opts = dict(width=750, height=350, fontsize={"legend": 8, "title": 10}, legend_position= "top_left")
            overtopping_levels = {"CN": 19.050, "CH": 18.840, "PSC": 25.5, "APC": 25.5} # unit: meters

            ## plot water level results
            cols = [col for col in df.columns if "WL_" in col]
            labels = [col.split("_")[1] for col in cols]
            plot_WL = hv.Overlay([hv.Curve(df[cols[i]], label=labels[i]).opts(tools=["hover"], ylabel="water level (m)") * 
                               hv.HLine(overtopping_levels[labels[i]]).opts(line_width=1.5, line_dash="dotted")
                               for i in range(len(cols))]).relabel("Results: water level").opts(**opts)

            ## plot discharge results
            cols = [col for col in df.columns if "Q_" in col]
            labels = [col.split("_")[1] for col in cols]                   
            plot_Q = hv.Overlay([hv.Curve(df[cols[i]], label=labels[i]).opts(tools=["hover"], ylabel="discharge (m3/s)") 
                                 for i in range(len(cols))]).relabel("Results: discharge").opts(**opts)
#             return (plot_WL + plot_Q).cols(1)
            return plot_WL, plot_Q

        except Exception as e:
            if sim.Mike1DException:
                return (pn.pane.Markdown(str(e), min_width=600, style={'color': 'red'}), None)
            else:
                raise e
                
    def view(self):
        obj1, obj2 = self.run_and_plot_results()
        if obj2:
            return (obj1 + obj2).cols(1)
        else:
            return pn.Column(pn.Spacer(height=20), obj1)
In [15]:
simulation = Mike1DSimulation_parameters(name="")
button = pn.widgets.Button(name='Run', button_type='primary', width=280)

def update(event):
    gui[1] = simulation.view()
       
button.param.watch(update, "clicks")
gui = pn.Row(pn.Column(simulation, pn.Spacer(height=50), button), simulation.view()).servable()

gui
# gui.show()
Out[15]:

Let's remember the inflows and the geosensor observations

In [16]:
plot_sensors_data.opts(width=500) + plot_inflows.opts(width=500)
Out[16]:

Calculate some basic statistics of the results

In [17]:
df_results = simulation.results
df_results.describe()
Out[17]:
WL_CN Q_CN WL_CH Q_CH WL_PSC Q_PSC WL_APC Q_APC
count 192.000000 192.000000 192.000000 192.000000 192.000000 192.000000 192.000000 192.000000
mean 18.309596 580.250987 18.199947 580.015704 18.155113 579.649339 17.653141 579.468992
std 0.547551 336.164841 0.497847 336.116437 0.478657 335.918510 0.112034 336.107876
min 17.680243 205.002643 17.642032 204.907977 17.626689 204.711055 17.528129 204.413269
25% 17.804565 279.353531 17.745624 279.978622 17.721412 279.729349 17.553057 279.053873
50% 18.272594 543.769740 18.155735 543.927904 18.106257 544.605283 17.642148 546.314709
75% 18.590993 732.706710 18.445257 732.418043 18.384517 733.388398 17.704425 733.353775
max 19.670065 1484.434576 19.470421 1485.963742 19.395857 1490.576550 17.953775 1481.509355